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ABSTRACT 

We analyze the mass distribution of cores formed in an isothermal, magne- 
tized, turbulent, and self-gravitating nearly critical molecular cloud model. Cores 
are identified at two density threshold levels. Our main results are that the pres- 
ence of self-gravity modifies the slopes of the core mass function (CMF) at the 
high mass end. At low thresholds, the slope is shallower than the one predicted by 
pure turbulent fragmentation. The shallowness of the slope is due to the effects 
of core coalescence and gas accretion. Most importantly, the slope of the CMF at 
the high mass end steepens when cores are selected at higher density thresholds, 
or alternatively, if the CMF is fitted with a lognormal function, the width of the 
log-normal distribution decreases with increasing threshold. This is due to the 
fact that gravity plays a more important role in denser structures selected at 
higher density threshold and leads to the conclusion that the role of gravity is 
essential in generating a CMF that bears more resemblance with the IMF when 
cores are selected with an increasing density threshold in the observations. 

Subject headings: ISM : clouds - ISM : globules - ISM : kinematics and dynamics 
- ISM : magnetic fields - turbulence - MHD 

1. INTRODUCTION 

It is now widely accepted that stars form in dense and gravitationally bound molecular 
cloud (MC) cores. Thus, understanding the properties of the core mass function (CMF) 
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is the cornerstone of any successful theory of the stellar initial mass function (IMF). The 
question whether there is a 1-to-l mapping between the CMF and the IMF remains an open 
one. The CMF is usually represented by dN/dM = M", where a takes different values 
in different mass ranges (for the IMF, a = —2.35 for M > 0.5 Mq, Salpctcr 1955). The 
determination of the CMF of cores of several MC regions have been obtained using a variety of 
wavelengths and techniques such as mm spectroscopy of various molecular lines, submm/mm 
continuum emission from dust, and stellar near infrared absorption by dust. Usually, no prior 
assumptions are made on the type of structures identified in these studies. They can be a 
mixture of prestellar cores, unbound objects, and eventually protostellar envelopes. For 
nearby MCs, where it is possible to identify dense cores (with average number densities 
n > cm~^; hereafter Cores), the CMF is usually described by a two component power 
law above and below a turnover point (with exponents and ai, respectively). Starting 
with the 1.2 mm study of the Ophiuchus starless dense cores by Motte et al. (1998), a 
number of submm/mm dust continuum observations have found ah and ai to be in the 
range [-2.1,-3] and [-1.3,-1.5], respectively (Motte et 2001; Testi & Sargent 2001; Johnstone 
et al. 2006; Beuther & Schilke 2004; Stanke et al. 2006; Reid & Wilson 2006; Johnstone 
& Bally 2006; Kirk et al. 2006; Enoch et al. 2006). These estimates bracket the values 
of the exponents of the IMF (Scalo 1998; Kroupa 2002; Chabrier 2003). Similar slopes are 
found using high density molecular tracers such as H^^CO"'" (Onishi et al. 2002; Ikeda et 
al. 2007) and the near infrared absorption technique (Alves et al. 2007). On the other 
hand, submm/mm observations of more distant clouds in which it is only possible to identify 
lower density structures (n ~ 10^ — 10^ cm~^; hereafter Clumps), exhibit values of in 
the range of [-1.4,-1.8] if the whole CMF is fitted with a single power law (Coppin et al. 
2000; Kerton et al. 2001; Tothill et al. 2002; Mookerjea et al. 2004; Massi et al. 2007; Li 
et al. 2007, Muiioz et al. 2007) or slightly steeper {ah ~ —2) if it can be fitted by a two 
component power law (Reid & Wilson 2005). Similar slopes are obtained from molecular line 
observations using low/intermediate density molecular line tracers (Stutzki & Giitsen 1990; 
Lada et al. 1991; Blitz 1993; Brand & Wouterloot 1995; Kramer et al. 1998; Heithausen et 
al. 1998; Kawamura et al. 1998; Wilson et al. 2003). Values of ah ~ [—1.6,-1.7] resemble 
the slopes of the stellar clusters mass function (Elmegreen 2000), which hints to the fact 
that these objects cannot be the direct progenitors of individual stars and must undergo 
additional fragmentation. 

The discrepancy in ah between the Cores and Clumps populations docs not seem to be 
merely an effect of the spatial resolution and hence of mass resolution. Additional evidence 
might be that, in Motte et al. (1998), the inclusion of the larger size cores in the CMF changes 
the values of from —2.5 to ~ —1.5. An explanation of the origin of the discrepancy would 
be that different physical processes dominate the evolution of these two populations of cores. 
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The close resemblance between the Cores CMF and the IMF suggests that they are more 
likely to be gravitationally bound (which is confirmed by mm line observations, Andre et 
al. 2007) and are the direct progenitors of stars or small stellar systems, whereas Clumps 
arc more likely to be dominated by turbulence. Note that the description of the CMF by 
power law functions is not the unique option. Reid el al. (2006) argued that a lognormal 
function is a better fit to the observed CMF in many regions and Goodwin et al. (2007) 
showed that the CMF of the submm cores obtained by Nutter & Ward- Thompson (2007) is 
well represented by a lognormal distribution (also Fig. 2 in Andre et al. 2008). 

Several physical processes arc believed to play a role in shaping the CMF: gravitational 
and thermal fragmentation (Larson 1985; Klessen 2001; Li et al. 2003), supersonic turbulence 
(Padoan & Nordlund 2002, PN02), accretion (Larson 1992; Bonnell et al. 2001; Bate & 
Bonnell 2005), core coalescence (Murray & Lin 1996; Dib et al. 2007a,2008), feedback 
(Adams & Fatuzzo 1996), and magnetic fields (Shu et al. 2004). A number of numerical 
studies of MCs that include some of these physical processes have derived a CMF that bears 
some resemblance with the CMF of Cores and the IMF (Gammie et al. 2003; Li et al. 
2004; Ballesteros-Paredes et al. 2006; Padoan et al. 2007) or with the CMF of Clumps 
(Hennebelle & Audit 2007, Heitsch et al. 2008). The origin of the differences between the 
Cores and Clumps CMFs and the eventual link between them has not been investigated 
so far. In this study, we identify cores in a high resolution 3D numerical simulation of a 
magnetized, turbulent, and self-gravitating MC model at different density thresholds. In 
Dib et al. (2007b), we have shown, using a detailed virial analysis that the inner parts 
of dense structures, identified at higher density thresholds, are more gravitationally bound 
than the same objects identified at lower thresholds. Thus, constructing the CMF at different 
thresholds is a direct test of the effects of gravity on the CMF. A striking example of how 
the density threshold can modify the inferred masses of cores is the case of the core L1512. 
Its derived mass from CO observations is ~ 6 Mq (Falgarone et al. 2001), whereas it was 
estimated from submm observations at ~ 1.4 Mq (Ward-Thompson et al. 1999). For the 
whole CMF, Kramer et al. (1998) derived a value of an = -1.72 ± 0.09 for the Orion B- 
S region from ^^CO(l-O) observations whereas Johnstone et al. (2006) derived a value of 
~ —2.5 from submm observations. Another example is the case of the NGC 7538 region, 
observed by Kramer et al. (1998) to have a CMF slope of an = -1.65±0.05 in the ^^CO(l-O) 
hne and ah — —1.79 ± 0.12 in the C^*O(l-0) hue which traces higher density material. 
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2. THE SIMULATIONS, CLUMP-FINDING ALGORITHM, AND 

ANALYSIS STRATEGY 

We performed a 3D numerical simulation of a magnetized, self-gravitating, and turbulent 
isothermal MC using the TVD code (Kim et al. 1999). The simulation is similar to the 
ones presented in Dib et al. (2007b) but with a higher numerical resolution of 512^ grid 
cells. In Dib et al. (2007b), we found that the population of cores formed in a moderately 
magnetically supercritical MC reproduces best some of the available observations such as 
the virial parameter-mass relation. Thus, we here only consider this cloud model for further 
analysis. The features of this simulation are: Turbulence is driven until it is fully developed 
(for 2 crossing timescales) before gravity is turned on. The Poisson equation is solved to 
account for the self-gravity of the gas using a standard Fourier algorithm. Turbulence is 
constantly driven in the simulation box following the algorithm of Stone et al. (1998). 
The kinetic energy input rate is adjusted such as to maintain a constant rms sonic Mach 
number Mg — 10. Kinetic energy is injected at large scales, in the wave number range 
k — 1 — 2. Periodic boundary conditions are used in the three directions. In physical units, 
the simulation has a linear size of 4 pc, an average number density of n = 500 cm~^, a 
temperature of 11.4 K, a sound speed of 0.2 km s~^, and an initial rms velocity of 2 km 
s~^. The Jeans number of the box is J^ox = 4 (number of Jeans masses in the box is 
Mbox/Mj^box = Jbox ~ 64, where Mhox = 1887 M©). The initial magnetic field strength 
in the box is Bq = 14.5 fiG. Correspondingly, the plasma beta and mass-to-magnetic fiux 
(normalized for the critical value for collapse M/0 = (47r^G)~^/^, Nakano & Nakamura 1978) 
values of the box are Pp^box — 0.1, and /ibox — 2.8. Cores are identified using a clump-finding 
algorithm that is based on a density threshold criterion and a friend-of-friend approach as 
described in Dib et al. (2007b). We restrict our selection of cores to epochs where the 
Truelove criterion (Truelove et al. 1997) is not violated in any of them. Thus, our derived 
CMF applies to starless prestellar cores. At a given threshold, the number of cores detected 
at one epoch is insufficient to sample appropriately the whole mass spectrum. It is therefore 
necessary to combine cores selected at several epochs in order to construct a CMF that 
is statistically significant. We combine the populations of cores identified at 35 timesteps. 
This approach resembles the situation in the observations where cores identified in different 
locations of a MC are most hkely in different stages of their evolution. Even for a composite 
multi-epoch spectrum, the statistics of cores identified with thresholds > 60 n remains poor. 
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3. CORE MASS FUNCTION 

Before discussing the CMF obtained from the simulation, it is worth checking what 
CMF is anticipated from the turbulent fragmentation formalism of PN02. Note that this 
formalism does not explain the CMF of Clumps. In the PN02 approach, ah is given by 
ah = [—3/(4-/3)] — 1, where /3 is the exponent of the turbulent velocity field power spectrum 
{Pvik) = k^^). In order to reproduce the CMF of Clumps, the formalism of PN02 requires 
/3 to be close to zero, which is not observed. Fig. [1] displays the time evolution of f3 in our 
simulation, and a few velocity spectra at selected epochs. The value of P is observed to 
fluctuate around ~ 1.9. According to PN02, this would result in a value of ah ~ —2.42 (at 
a given epoch and for a multi-epoch built CMF). Fig. [2] (left) displays the CMF for cores 
identified at the density thresholds of 20 and 50 n. Both samples use the same number of 
timesteps. Several features of the CMFs can be noted. If the high-mass end of the spectra is 
to be assimilated to a power law, fits to the CMF in this regime yield a steeper slope for the 
high density threshold sample than for the one with the lower threshold. The slopes for the 
20 n and 50 n CMFs are —1.55 ± 0.15 and —1.80 ± 0.2 in the mass ranges [1,113] Mq and 
[2.8,44] M0 respectively. However, the CMFs are also well fitted over their entire mass range 
by a lognormal distribution whose width is smaller at higher threshold. The width of the 
lognormal fit of the CMF is 1.63 ±0.07 and 1.12 ±0.09 at the 20 and 50 n levels, respectively 
(the width of the IMF is 0.57, Chabrier 2003). A two-sided Kolmogorov-Smirnov (KS) test 
indicates that the 20 n and 50 n CMFs have a probability of only 2 percent of being drawn 
from the same parent distribution, a result significant at the 2.3 sigma level. For reference, we 
also calculated the CMF slope of 20 n cores in a non self-gravitating model (Fig. [2l right). 
The slope of the CMF in the high mass range is larger for the non self-gravitating cloud 
model than for the self-gravitating one {ah = —1.73 ± 0.07 and resembles the one predicted 
and obtained by Hennebelle & Audit (2007) for turbulence dominated, non self-gravitating 
clumps (in this non self-gravitating case, statistics are not enough to build the 50 n CMF). 

The slopes of the CMF at high mass end are shallower than those predicted by PN02 
and steepen with increasing threshold. We interpret the shallower than the IMF slopes 
and the steepening with increasing threshold to the role played by gravity in modifying 
the CMF generated by turbulence. The steepening of the slope with increasing threshold 
is due to the increasing importance of gravity at higher thresholds. Fig. [3] displays the 
ratio of the gravitational term in the virial theorem [W = — jy pri{d(j)/ dri)dV , where (f) is 
the gravitational potential) to the other terms that appear in the virial equation. The latter 
terms are condensed in the quantity Qvt = 2 {Eth + Ek — Tk — Tth)+Emag + 'Tmag, where Eth = 
I fy PdV, Ek = \ ^y pivfdV, and Emag = /y ^^^V are the volume thermal, kinetic, and 
magnetic energies, respectively, with P being the thermal pressure, and Tth = ^ ''^iPnidS, 
Tk = \ jgTipViVjhjdS, and r^ag = ^i^ij^c^'S' are the thermal, kinetic, and magnetic surface 
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energies, respectively, where Tij is the Maxwell stress tensor. Fig. [3] indicates that gravity is, 
on average, more important in the virial balance for the cores defined at the 50 n threshold 
which are a mixed population of gravitationally bound and unbound objects than for the 20 n 
cores which are essentially unbound objects. An extrapolation to higher thresholds of the 
order of 500 n = 2.5 x 10^ cm~^ which are number densities similar to those of gravitationally 
bound cores observed in the submm/mm or in dense molecular tracers would yield slopes 
in the range [-2.3,-2.5], thus offering a natural explanation for the similarity of the observed 
slopes between the CMF of dense cores and of the IMF. The reason why the CMF of the 
20 n cores in the self-gravitating regime is shallower than that of the same population in 
the non self-gravitating case is most likely due to the long range effects of gravity. Fig. H] 
displays the time evolution, for the 20 n cores, of the fraction of the mass present in the 
cores, number of cores, and average and median masses of the cores. The increasing mass 
fraction and average mass of the cores at roughly equal number of cores even before any of 
them proceeds to collapse are a clear indication of mass growth affecting the cores due to 
the combined effect of gas accretion and core coalescence (see Fig. 16 in Dib et al. 2007b). 
On the other hand, the constancy of the number of cores and their median mass means 
that newborn cores are generated with the smallest masses. Although it is not possible to 
assess from Fig. H] the contribution of each of these two processes, this mass growth which 
constitutes a deviation from the purely turbulent regime is responsible for the flattening in 
the 20 n CMF in the self-gravitating case as compared to the non self-gravitating one. 



4. CONCLUSIONS 

In this work, we have constructed the mass function of cores in a 3D numerical simulation 
of a magnetized, turbulent, and self-gravitating molecular cloud. Cores are identified at 
the density threshold levels of 10^ and 2.5 x 10^ cm~^. In agreement with observations, 
the slope of the core mass function is found to be steeper for cores identified at the higher 
threshold level. The steepening of the slope with increasing threshold is due to the increasing 
importance of gravity in the virial balance of the cores at higher densities and a sign of the role 
played by gravity in the gravo-turbulent scenario of star formation. On the other hand, at 
lower thresholds, the slope is shallower than what is found for the purely turbulent case. This 
is due to the effects of mass growth that affect the cores. This mass growth is essentially due 
in the initial phase to the effect of core coalescence which generates more massive structures 
particularly for cores defined at lower thresholds since they have larger cross sections. 
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Fig. 1. — Time evolution of the exponent of the turbulent velocity field power spectrum, /? 
(top) after gravity is turned on, and a sample of velocity power spectra at selected epochs 
(bottom). Each data point in the top figure is an average over five neighboring timesteps. 
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Logio(M/M3„,) Logio(M/M3„,) 

Fig. 2. — Left: The CMF of cores identified in the simulation with gravity at the thresholds 
of 20 n (thick full line) and 50 n (thin dash-dotted line). Overplotted are lognormal fits (thin 
line) and power law fits to the high mass end (dashed line). Right: The CMF of cores at 
the threshold of 20 n for the models with gravity (thick full line) and without gravity (thin 
full fine) . Power law fits to the high- mass end are overplotted (dahsed fines) . 
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Log(IWl/l0vTl) 

Fig. 3. — Normalized probability distribution function of the ratio of the gravitational term 
W to both volume and surface energy terms that appear in the virial equation (kinetic, 
thermal, and magnetic; full hne is for the 20 n cores and the dashed line is for the 50 n 
cores). The latter terms are condensed in the quantity Qvt (see text for details). The 
larger |H^|/|0vt| is, the stronger the influence of gravity. A KS test of the (non-normalized) 
distributions shows that they have a probability of 3 percent of being identical. 
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Fig. 4. — Time evolution of: The fraction of the mass present in the cores (top, full line), 
number of cores (top, dashed line), average mass of the cores (bottom, full line), and median 
mass of the cores (bottom, dashed line). Time is shown in Myrs (bottom horizontal axis) 
and in units of the crossing time tc— 2 Myr (top horizontal axis) . 



